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It is known that ultrasound techniques yield non-intrusive measurements of hydrodynamic flows. 
For example, the study of the echoes produced by a large number of particles insonified by pulsed 
wavetrains has led to a now standard velocimetry technique. In this paper, we propose to extend 
the method to the continuous tracking of one single particle embedded in a complex flow. This gives 
a Lagrangian measurement of the fluid motion, which is of importance in mixing and turbulence 
(~| ■ studies. The method relies on the ability to resolve in time the Doppler shift of the sound scattered 

by the continuously insonified particle. 
I For this signal processing problem two classes of approaches are used: time-frequency analysis and 

I , parametric high resolution methods. In the first class we consider the spectrogram and reassigned 

^ ' spectrogram, and we apply it to detect the motion of a small bead settling in a fluid at rest. In more 

' non-stationary turbulent flows where methods in the second class are more robust, we have adapted 

^ , an Approximated Maximum Likelihood technique coupled with a generalized Kalman filter. 

CJ 

• ^ ■ PACS numbers: 43.30.Es, 43.60.-c, 47.80.-l-v, 43.60.1 

Qj;_ I. INTRODUCTION 

' In several areas of fluid dynamics research, it is desirable to study the motion individual fluid particles in a flow, 
^ ' i.e. the Lagrangian dynamics of the flow. The properties of this motion governs the physics of mixing, the behavior of 
' binary flows and the EuleriAn complexity Qfekaptic and turbulent flows. Lagrangian studies are possible in numerical 
25 . experiments where chaoticEl and turbulenta'ElQB flows have been studied. For turbulence, the numerical studies are 
limited to small Reynolds number flows whose evolution is only followed during a few large-eddy turnover times. In 
addition only the small scales properties of homogeneous turbulence are captured; the influence of inhomogcneitics 
(such as large scale coherent structures) are not taken into account. It does not seem possible at the moment to extend 
high resolution turbulent DNS computations to long periods of time or to high Reynolds number flows. Experimental 
studies are thus needed. They differ from the numerical studies because one cannot tag and follow individual fluid 
particles; most techniques aim at recording the motion of solid particles carried by the flow motion. The degree 
' ^ ' fidelity with which solid particles can act as Lagrangian tracers is an open problem; it depends on the size and 
J>-,' density of the particle. While the interaction between the particle atyiits wake can be important for large particles or 
. particles with a large density difference with the surrounding fluidOulj, it is generally admitted that density matched 
particles with a size smaller that the Kolmogorov length follow the fluid. Measurements of small particle motion have 
been made, using optical techniques that follow individual particle motion over short times/distanceaail3. We propose 
here an acoustic technique that can resolve an individual particle motion over long periods of time (compared with 
the characteristic time of flow forcing) . 
^ . The principle of the technique is to monitor the Doppler shift of the sound scattered by a particle which is con- 
tinuously insonified. This is an extension of the pulsed Doppler method tlaat has been developed to measure velocity 
profiles and that has many applications in fiuid mechanics and medicinetj. The main advantage of the continuous 
insonification is to improve the time resolution of the measurement, although it is limited to the tracking of a very 
small number of particles (the tests reported here are made with only one particle in the flow). The measurement 
relies on the ability to track a Doppler frequency and its variation in time. For this signal processing problem two 
classes of approaches have been developed: (i) time-frequency analysis and (ii) high resolution parametric spectral 
analysis. Time frequency methods mainly rely upon the quadratic Wigner-Ville transform, or smoothed versions of it. 
Numerous studies ancL^apers have-recently been published, in which the theoretical issues are presented (see e.g. the 
textbooks by Flandrinll3 or Cohenll3) . These non parametric techniques are convenient and well-suited for weakly non- 
stationary signals with a good signal-to-noise ratio (SNR). However, time frequency representations present numerous 
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FIG. 1; Principle of measurement. A large 3D measurement zone is achieved by using a transducer size of a few wavelengths. 



drawbacks when it comes to extract trajectory information. Their quadratic nature give rise to numerous spurious 
interference terms that require post processing. For signals with a faster frequency modulation and a low SNR, we 
show here that an optimized parametric approach is a better choice. Parametric high resolution spectral analysis 
methods take advantage of an a priori knowledge of the spectral content of the recorded signal, namely the emitted 
signal frequency plus one or many doppler-shifted echoes. Furthermore, a time-recursive frame for the estimation of 
the Doppler shift is proposed here, where the evolution of the frequency is taken into account in the algorithm. 

The two methods are tested in two experiments, in which the acoustic signals have different time scales and noise 
levels. The first experiment is a study of the transient acceleration of a heavy sphere settling under gravity in a fluid 
at rest. In this case the characteristic time scale of velocity variations is slow (r ^ 50 ms) and the signal to noise ratio 
is fair (about 20 dB); we show that a technique of reassignment of the spectrogram gives good results. The second 
experiments deals with the motion of a neutrally buoyant sphere embedded in a turbulent flow. In this case, velocity 
variations occur over times of about 1 ms and the signal to noise ratio is low (less than 6 dB). We show that the AML 
parametric method yields very good results in that situation. 

The paper is organized as follows: in section II we present the acoustic technique and measurement procedure. In 
section [II we describe the signal processing techniques, with a particular emphasis on the AML method which has 
been developed and optimized to this particle tracking problem. Examples of applications to measurements in real 
flows are given in section IV. 



II. ACOUSTICAL SET-UP 



A. Principle of the measurement 



In the experimental technique proposed here, a particle is continuously insonifled. It scatters a sound wave whose 
frequency is shifted from the incoming sound frequency due to the Doppler effect. This Doppler shift is directly 
related to the particle velocity v^: 

Aw = q-Vp , (1) 

where q is the scattering wavevector (the difference between the incident and scattered wavevectors q = kgcat ~ ^inc) 
and uj is the wave pulsation. 

We choose a backscattering geometry (see figure |l|) so that q = — 2kinc and the frequency shift becomes 

ALj{t) ^ -2^LOo , (2) 
c 

where c is the speed of sound, ujo is the incident pulsation, and v(t) is the component of the velocity on the incident 
direction at time t. We continuously insonify the moving particle and record the scattered sound. If need be, the 
particle position can be obtained by numerical integration of the velocity signal. 
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FIG. 2: Data from a steel bead (diameter 1 mm) settling in water at rest, (a) Typical time series; (b)power spectral density of 
the inset figure. On the x-axis, zero corresponds to the emission frequency. 

B. Transducers characteristics and acquisition 

We use a Vermon array of ultrasonic transducers made of individual elements of size 2x2 mm each, separated by 
100 /im. Their resonant frequency is about 3.2 MHz and their bandwidth at -3 dB is 1.5 MHz. Sound emission is 
set at 3 MHz or 3.5 MHz; experiments are performed in water so that the wavelength is A=0.50 mm or 0.43 mm. 
The corresponding emission cone for each d = 2 mm square element is 29° at 3 MHZ and 24° at 3.5 MHz. In our 
measurements, the particle to transducer distance lies between 5 cm and 40 cm, so that measurements are made in 
the far field {(P /\ > 10 mm). Given maximum flow and particle velocities of the order of 1.5 m.s~^, we expect a 
maximum sound frequency shift of the order of 5 kHz or 6 kHz, depending on whether the emission is at 3 MHz or 
3.5 MHz. This yields a frequency modulation rate of at most 0.25%. One element of the transducer array is used for 
continuous sound emission and another for scattered sound detection. As the operation is continuous (as opposed to 
pulsed) and the elements are located close to one-another, we observe a coupling between the emitter and the receiver 
of the order of 60 dB (this is due both to electromagnetic and acoustic surface waves cross-talk). 

The sound scattered by the moving particle is detected by a piezoelectric transducer. Upon connection to a 50 
impedance, it yields an electrical signal of about 2 to 30 /iV. In comparison, the noise is 1 /^V and the electromagnetic 
coupling with the emitter is 8 mV. Hence the signal to noise ratio is between dB and 30 dB. The transducer output 
is sampled at 10 MHz over a 21 bit dynamical range (input range 31.25 mV) and numerically heterodyned at the 
emitting frequency. Then it is decimated at the flnal sampling frequency of 19531 Hz. The acquisition device is a 
HP-el430A VXI digitizer. 



C. Scattering by an elastic sphere 

The study of sound scattered by a flxed solid sphere is a classic but continuing area of study and dif 
arise in the interpretation of observed phenomena especially when trying to deal with elasticity and absorption 
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Complex behaviour is observed linked with resonances of Rayleigh waves at the surface of the sphere. As a consequence 
the scattered pressure distribution varies both in directivity and amplitude. A generic expression for the far field 
pressure is the following : 

Pscat{r,6)^p.^,'^i^^e^^^ , (3) 
Zr 

where r is the distance from the center of the sphere, a its radius, pinc the incident pressure on the sphere, k the 
incident wavenumber in the fluid, the scattering angle and / is a form function which depends on the physical 
properties of the solid medium. Under very general assumptions, / can be developed as a series of partial waves: 

/(^«'^) ^ l|^E/2n+ l)|iMp„(cos0) , (4) 

where P„ is the a Legendre polynontial, i?„ and D„ are determinants of matrices composed of spherical Bessel and 
Hankel function and their derivativeao. Physically, / represents the sum of the specular echo and of interferences due 
to the radiation by Rayleigh waveaijO. As a result, / is a strongly varying function, particularly for high values of 
ka. In our experiments we used spheres of different material (polypropylene PP, steel, tungsten carbide, glass) with 
corresponding ka between 7 and 15. The flow acts on the sphere motion, thus causing its acceleration and, eventually, 
its rotation. These effects may change the radiation diagram: first there is Doppler shift for the sound received by 
the sphere, and, perhaps more importantly, the sphere rotation may change the Rayleigh emission. For these reasons, 
the evolution of the amplitude of the scattered sound during the particle motion is quite complex. However, the 
observed amplitude modulation (see figures ^ and ||) varies slowly enough to allow a correct estimate of the frequency 
modulation of the scattered sound. 



III. SIGNAL PROCESSING 



Numerous spectral estimation techniques are based on the ideas behind Fourier analysis of linear time invariant (LTI) 
differential equations. These techniques may be divided into (i) non-parametric techniques where the basis functions 
are implicitly the harmonically related complex exponentials of Fourier analysis and (ii) parametric techniques whose 
task is the estimation of the parameters of a (sub)set of complex exponentials. The spectrogram and the reassigned 
spectrogram belong the former category, whereas the maximum likelihood and its approximate form belong to the 
latter. 



A. Time-Frequency analysis 

The most common time frequency distribution (TFD), the spectrogram, involves a moving time window. This 
window attempts to capture a portion of the signal which is sufficiently restricted in time so that stationarity and 
LTI assumptions are approximately met. TOjeyercome the inherently poor localization in the time-f rc a u ency plane, a 
method has been proposed by Gendrin et al.LD, and extended more recently by Auger and Flandrinll3li3. The idea is 
to locally reassign the energy distribution to the local center of gravity of the Fourier transform. Despite its ability to 
exhibit clear and well localized trajectories in the time-frequency plane, this technique requires an additional image 
processing step to extract the TF trajectory. For rapidly fluctuating frequency modulations and/or low SNR spurious 
clusters appear which makes this extraction difficult. The parametric method presented below is more robust. 



B. AML spectral estimation 

This approach is largely based upon maximum likelihood spectral estimation (see e.g. KayEi). The fundamentals 
are briefly rejcalled, as they serve as a basis for the approximate likelihood scheme, originally developed by Clergeot 
and TressensEEI. This work is extended here within a recursive estimation frame, thus allowing to track the variations 
of the Doppler frequency shift induced by fast velocity changes of a scattering sphere imbedded in a turbulent flow. 
Michel and Clergeot have developed a similar approach for non stationary spectral analysis in an array processing 
framcEllEj. 
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1. Introduction 



In this section, we address the problem of estimating the frequencies /i , . . . , /m of M harmonic signals embedded 
in noise, from a small number of samples 

M 

x{t) = ^ a™(i) exp(j(27r/™i + + n{t) . (5) 

m—l 

As the number of sampling points that are supposed to be available is low, classical Fourier based approaches fail 
to provide good results. We focus our attention on parametric approach, where an a priori knowledge about the 
structure of the signal is taken into account to improve the analysis. 
The following assumptions are made: 

• The time series is regularly sampled with time period Tg, so as to insure ^ > ^"'2" where fmax stands for the 
bandwidth of the anti-aliasing filter used in the recording process. For convenience, T, will be set to = 1 and 
the term frequency will refer to normalized frequency (i.e. the actual frequency, divided hy Fg — . 

• The amplitudes ai{t) are deterministic but unknown. 

• The noise is a complex white gaussian circular and iid (independent increment identically distributed) with 
(unknown) variance u^; the distribution function of a fc-dimensional vector N defined by 

N{t) = [n{t),n{t + 1), n{t + 2), . . . ,n{t + {K - l))f (6) 

reads 

Furthermore, the noise and the signal are independent. 

• The term observation refers to a set of Q iC-dimensional vectors constructed from the sampled time series, 
according to 

X(t,) = [x{t,),x{t, + l),x{t, + 2), . . .,x{t, + {K- l))f J = 1, . . . , Q (8) 

• The frequency / is supposed to remain constant during an observation. 

Under the assumption that the noise process is iid, and that the observed vectors are corrupted by independent reali- 
sations of the noise process, the likelihood of the observation is given by the product of the likelihood of each vector. Let 
V be the set of searched parameters {V contains a^, F = {/i, . . . , f]\f} and the A = [ai exp{j(j)i), . . . , om eyip{j(j)M]'^), 
the loglikelihood of an observation is simply given by 

1 

C{V) = -KQ \og{2na') - ^ E " S(F)A(q)p , (9) 

^ 9=1 

where 

/ 1 exp(27r/i) . . . exp(27r(ii: - l)/i 



S(F) = 



(10) 



V 1 exp(2^/M . . . exp(27r(i^ - 1)/, 



M 



According to the maximum likelihood principle, the set V of parameters must be chosen in order to maximize 
expression (0). 
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2. Reduced expression 

Minimizing jointly for all the parameters is usually untractable. Most authors propose a separate maximisation 
for each of the parameters. For our application, the spectral components (i.e. A and F) are the relevant variables. 
We first maximize with respect to A and derive an expression for the optimal F; is estimated independently. 

The value of vector A which minimizes the norm \X{q) — S{F)A{q)\'^ is easily obtained : 

Aiq) = (S+S)-iS+(F)X(g) . (11) 

Note that the "signal only" vector Y — X — N appears to be the orthogonal projection of X on the signal subspace 
spanned by the row vectors of S: 

Y = S{F).A{q) = S((S+S)-iS+(i?)X = Us{F)X , (12) 

where ns(-F) stands for the parametric projector on the signal subspace0. Let n„(i?) = I — ns(F) be the noise 
subspace, I is the identity matrix. By substituting ( |l^ and using the definition of n„(_F) in the expression of the 
log-hkelihood (|), one gets the following simplified expression to minimise 



1 

L{F)^-Y,\n,,{F)X\\ (13) 

q=l 

Using the properties of the trace operator (hereafter denoted Tr) and those of the projection matrix n„(i^), the 
maximum likelihood estimation of F takes the common form : minimize 

L(F) = -^Tr [n„(F)RJ , (14) 

where R^; is an estimate of the correlation matrix Rj; of the vector process X{q). Minimizing L(F) in (|lj) leads to the 
exact value Fml which has the maximum likelihood. It is important here to emphasize the following : if the vectors 
X are obtained by time-shift over the recorded time series, an observation runs over K + Q — 1 samples, i.e. the actual 
duration of one observation is Tots = {Q + K — 2)Ts. In this case, the observed vectors may not be considered as being 
corrupted by independent realisations of the noise process, as some 'time integration' is performed iif-the estimation_af 
Tlx- The consequences and interest of such smoothing have been studied by Clergeot and TressenscB, and OuamriEj, 
in the frame of array processing (in this context, 'time integration' becomes 'spatial smoothing'). In the remainder 
of this paper, the development are based on equation (p^, no matter how R^; is estimated; see appendix for the 
practical implementation. 



3. Approximate Max-likelihood 



Equation ( [l4| ) is still too complicated to be solved analytically in a simple way. A minimization can be easily 
performed if L{F) has a quadratic dependance in Sc3. Let Ry be the correlation matrix of the signal vectors Y{q), 
the assumption that signal and noise are independent allow to establish the following equalities 



Rx — Ry -\- (J 1. 

Ry = SPS+ 



T>^£[AA+] , 

where £ stands for the mathematical expectation. Substituting in equation (p^) leads to: 



L(F) 



rTr 



^„(F)SPS^ 



Clergeot and TressensB propose a second order approximation of L{F): 



Laml(F) = ^Tr 



n„S(F)PS+(F) 



(15) 
(16) 
(17) 



(18) 



(19) 
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in which ILn is estimated by computing the projector spanned by the {K — N) smallest eigenvalues of the estimated 
covariance matrix R^, . They prove that this approach leads to more reliable estimates of F at low signal to noise ratio 
(SNR), and that the minimization of Laml is asymptotically efficient. In practice, the following set of equations is 
used 

= ^^^'■(^nR.) , (20) 

n,(F) = S(F)(S+(F).S(F))-i.S+(F) , (21) 
S(F).P.S+(F) = n,(F)(R, - a2l).n,(F) . (22) 

The approximately quadratic dependence of Laml in S(_F), allows a fast convergence of the minimization algorithm 
by using a simple Newton-Gauss algorithm: 

F(k + 1) = F{k) - H-i.grld(i^ML)|j?^j?(fe) , (23) 

where k stands for the iteration step in the minimization process, grad and H are the gradient and hessian respectively 
(see expressions in the appendix). 



4.. Combining new measurements and estimates 

In this section, it is assumed that new measurements do not allow by itself the derivation of a good estimate. The 
variance of such an estimate varies as y~ — + Q ~ l)~^j whereas integrating new measurements to this estimate 
allow to derive a better estimation. ^ 

Let F{t) be an estimate of F at time t, and Af{F{t),T{t)) its density, assumed to be normal with variance rit)M. If 
a linear evolution model is known for F{t), one has 

F(t + 1) = MF{t) + e{t) , (24) 
Pt+i\t{F) = J^{MP{t), Mr(i)M+ + R,) , (25) 

where M is the evolution matrix; e is a perturbation term, which is statistically independent from F, and , Rg is its 
covariance matrix. Pt+i\t is the probability density function that can be derived for time i + 1, if the observations are 
made until time t only. As sucL-an evolution equation is usually unknown, M will be set to the identity matrix in 
the rest of the paper (see MichelE3) for a detailed discussion). Applying the Bayes rule over conditional probabilities 
gives : 

. Pt+^t{F).pt+i{X\F) 

Pt+i\t+i[F) = . (26) 

Pt+i{X) 

Noting that log(pt+i{X\F)) is the loglikelihood function for which a reduced expression has been derived in the 
previous section, one gets after all reductions and identifications the simple following expressions 

Fit + l\t)=F{t) , (27) 
r(i + l|t) = r(i) + Re , (28) 

r(t + i)-i = H + r(t + i|f)-i , (29) 

F{t +l\t+l) = F{t + l)= F{t + l\t)- T{t + l)-i.grad , (30) 

where it can be shown that the gradient function has the same expression as in the previous section. Rg is an unknown 
matrix which will be practically set to w^I, where v'^ will be tuned in order to allow the algorithm to take slight changes 
in F into account. Furthermore, it is interesting that the set of expression above expresses a generalized Kalman 
filter for estimating F (in the sense that it relies upon second order expansion of the loglikelihood functions). The 
statisticaJ-convergence^roperties and numerical efficiency of these approaches are described in the work of Michel & 
ClergeotEl and MicheE3. 
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FIG. 3: Experimental setup in the case of the settling sphere. 

IV. EXPERIMENTAL RESULTS 

We first describe the simple case of a particle settling in a fluid at rest. It is well adapted to the reassigned 
spectrogram method because the acoustic signal has a good SNR and a slow frequency modulation. We show that it 
allows to extract the subtle interaction between the falling particle and its wake. We then study the more complicated 
case of the motion of a particle embedded in a turbulent flow, where the dynamics of motion is much faster and the 
SNR is poor. We show that the AML method is well suited. 

A. The settling sphere 

1. Motivation and experimental setup 

When a particle is released in a fluid at rest, its developing motion creates a wake. The particle velocity is then 
set by the balance between buoyancy forces and drag, and additional subtle effects: first, 'added mass' corrections 
because the particles 'pushes' the fluid, and seco]ad,|-a 'history' force because the wake reacts back on the particle. 
Formally, one can write the equation of motion asoBo: 

{mp + f)^^ = {jUp- mf)g- ^na^Pf \\vp\\vpCD{Re) + FhistoTy ■ (31) 

where rup is the particle mass, is the mass of a fluid particle of the same size, Vp is the particle velocity, g is the 
acceleration of gravity, a is the sphere radius, pf is the fluid density, cd is the static empiric drag coefficient, Re is 
the Reynolds number Re = and Fhistory is the so-called history force. In this expression, the drag coefficient is 
usually obtained from measurement of the forces acting on a body at rest in an hydrodynamic tunnel. The history 
term, however, is largely unknown. Analytic expression can only be derived in the limit of small Reynolds numbers 
(less than 10) and cannot be applied for real flow configurations (e.g. multiphase flows) where Re ^ 1. 

We perform measurements of the motion of a settling sphere, with the aim of evaluating the influence of the history 
forces. We use a water tank of size 1.1 mxO.75 m and depth 0.65 m, filled with water at rest (figure The bead is 
held by a pair of tweezers, five centimeters below the transducers. It is released a time t=0 without initial velocity 
and its trajectory is about 50 cm long. The data acquisition is started before the bead is released in order to capture 
the onset of motion. 



2. Results 

Let us use as a first example, the fall of steel bead, 0.8 mm in diameter. The Doppler shift during the bead motion 
is detected using the spectrogram representation and a subsequent reassignment scheme. The simple spectrogram 
and reassigned version are shown in figure IJ. The reassignment technique drastically improves the localization of the 
energy in the time-frequency plane. In this case, the image processing step computes Vp{t) as the line of maxima. 
The precision of the overall measurement depends on two factors : first on the intrinsic precision of the reassignment 
method and second on the dispersion of the measurements (the reproducibility of the bead motion over several 
experiments) . The intrinsic precision of the reassignment method has been empirically studied using synthetic signals 
modelling the particle dynamics plus a noise that mimics the experimental data. We observed that for our choice of 
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FIG. 4: (a) Spectrogram of the backscattered sound, after heterodyne detection. (6) Reassigned spectrogram. In each figure 
the inset shows a norrrpiized cross-section of the spectrogram. The algorithm is that of the tf rrsp function of the MATLAB 
time-frequency toolboxO. To get rid of the spectral components at zero frequency due to the coupling between transducers 
and at small frequencies around zero due to slow motion of the water surface, we use a high pass fifth order Butterworth filter 
of cut-off frequency 25 Hz (corresponding to a velocity of 5 mm/s). Data of 0.8 mm steel bead settling in water at rest. 



parameters (a time- frequency picture with 256x256 pixels) the rms precision is about one half pixel both in time and 
frequency directions. The method thus allows a precise analysis of the dynamics of the fall; we describe below two 
sets of experiment that illustrate the potential of the reassignment technique. 

First, we show in figure ^ the velocity of a 1 mm steel bead (average over ten falls) together with two numerical 
simulations based on equation ^l], first without the memory force and second with the expression of the memory 
force derived at low Reynolds numbers (called the Stokes memory term, as in Maxey & RileyH). The precision of 
the detection technique is sufficient for the measured profile to be compared to the simulated curves and to draw 
physical conclusions about the hydrodynamical forces. At early times, the trajectory is close to the simulation_with 
memory force. This is due to the diffusion away from the bead surface of the vorticity generated at the boundaryElBcj. 
However, as the instantaneous Reynolds number increases, the curve deviates from this simple regime: vorticity is 
advected into the wake. Memory is progressively lost and the sphere reaches a terminal velocity in a finite time as 
does the simulation without memory. 

The measurement and signal processing techniques are then tested on a more non-stationary motion, as in the case 
of a bead whose density is close to that of the fluid. In this situation a stronger interaction is expected between the 
particle motion and the development of its wake. Formally, this traces back to differences in the effective inertial mass 
and buoyancy mass of the particle - see equation (p^) . In figure ^, we show the velocity variation for a light glass 
sphere (density 2.48) compared to a tungsten bead (density 14.8). We observe that the velocity of the glass oscillates 
before reaching a constant terminal value whereas the other particle has a regular acceleration. In the case of light 



10 




FIG. 5: Velocity measurement of a steel bead of diameter 1 mm (solid line), compared to numerical simulations without memory 
force (dashed) and with Stokes memory (dash-dotted). The inset shows an enlargement near the onset of motion. The Reynolds 
number, based on the limit velocity is 430. The sphere velocity profile results from averaging n = 10 successive experiments. 
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FIG. 6: Fall of a tungsten carbide sphere D=l mm (dashed) compared to a glass bead D=2 mm (solid), at Re ~ 400. The 
velocity is non-dimensionalized by the limit velocity. Curves are not averaged over several experiments. 



beads the hydrodynamic forces may be large enough to overcome the gravity and change the sign of the acceleration 
. This is linked with the non-stationarity of the wake, as vortex shedding is known to occur for Reynolds number 
above critical {Rec ~ 250). 



B. Turbulent flow : Lagrangian velocity measurement 



1. Experimental set-up 



The turbulent flow is generated in a von Karman geometry : the water is set into motion by two coaxial counter 
rotating disk in a cylindrical tank (figure |^). The Reynolds number Re = ^ (where v = 0.8910"^ m^s~^ is the 
water kinematic viscosity) is equal to 10^. To prevent cavitation in the flow, we boil the water before filling the tank 
by lowering the pressure with a vacuum pump and during the experiment the pressure is increased to two bars. For 
the acoustic measurement, we use the same array of transducers as in the previous experiments, at emitting frequency 
3 MHz. The cylinder and the surface of the disks are covered by 3 cm of Ciba Ureol 5073A and 6414B. Its density 
is 1.1 and the sound velocity is 1460 m.s~^ so that its acoustic impedance is close to that of the water, reducing 
drastically the reflections at the interface water/ureol compared to water/steel. The attenuation at 2.5 MHz is about 
6 dB per cm. With a 3 cm layer and after the reflection on steel the total absorption is about 36 dB. The total 



11 




FIG. 7: Experimental setup. The inner radius of the cylinder is 10 cm (disks radius R—9.5 cm) and the distance between the 
disks is 18 cm. The disks are driven by two 1 kW motors at a constant rotation frequency of /=18 Hz. The transducers are 
placed 18 cm off axis, in order to increase the volume of the measurement region. 



reflection at the interfaces is reduced by a factor 60. The particle is a polypropylene (PP) sphere of radius 1 mm and 
density 0.9. 

2. Results 

We show in figure || the time series when one particle is in the ultrasonic beam and the corresponding spectrogram 
and reassigned spectrogram. The signal to noise ratio is very poor, typically less than 6dB (to give an idea, in 
figure ^(a) the bead enters the ultrasonic beam at i ^ 20 ms) . One can also see some events localized in the time 
frequency plane that may be considered as noise and that may have several origin (noise of the motors, external 
electromagnetic noise, cavitation in the flow. . . ). Altogether, the time-frequency pictures show the trajectory of the 
particle but the low SNR prevents it from being easily extracted. In particular, the trajectory in the reassigned picture 
becomes quite lacunar and extracting it would require sophisticated (and CPU greedy) image processing techniques. 

The result of the AML algorithm is plotted in figure ^. The extracted frequency modulation is of course within the 
estimation in the spectrogram as in figure H(b), but one observes that fine variations in the velocity of the bead are 
now detected. The algorithm provides also an estimate of the amplitude of the source (figure ^(c)). It can be seen 
that there is a strong amplitude modulation and that the SNILjs at most 6 dB and may become less than dB. 

As the hessian is related to the Fisher information matrixes, its inverse square root is linked with the variance 
of the estimation: a large value of the hessian indicates an accurate estimation of the modulation frequency and, 
hence, of the bead velocity. The inverse square root of the hessian is plotted in figure ||(b): very large values are 
calculated in the absence of a bead in the measurement volume at the beginning and end of the time series (as a 
signature of the mismatch between the model which is composed of one source at least and the reality: no source) . 
Local lower values (typically less than 0.1) are observed when the variance on the estimation is small. Spurious effects 
are generated when the frequency modulation approaches zero as the hessian also becomes very small because of the 
filtering operation made in order to get rid of the coupling part of the signal. Finally, one observes that the hessian 
decreases as the signal to noise ratio increases (see at time 0.55 s). 

V. CONCLUDING REMARKS 

As can be seen in the previous section, both methods, time-frequency analysis and parametric spectral analysis are 
suited for extracting the time- varying frequency modulation due to a Doppler effect. The domain of application of 
each method depends on the degree of non-stationarity and on the SNR. 

For high SNR and weakly non-stationary signals, the time-frequency approach yields very good results. One 
drawback is the need of a second processing stage to extract the trajectory from the time-frequency picture. This 
stage may become increasingly difficult if there is more than one spectral component or if the SNR degrades. In 
both cases the quadratic nature of the algorithm produces interference patterns in the image: spurious clusters and 
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FIG. 8: Sound scattered by a 2 mm diameter PP bead in a turbulent flow at Re = 10^. (a) Typical time series; (b) and (c) 
corresponding spectrogram and reassigned spectrogram. 
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FIG. 9: Velocity measurement for the motion of a 2 mm diameter PP bead in a turbulent flow at Re = 10®. Output of the 
AML algorithm: (a) velocity (b) corresponding inverse square root of the hessian, (c) amplitude of the source (the rms value 
of the noise is 0.9//V). AML algorithm parameters: M=l, K=7, Q=13, v^=10"^. 
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a lacunar trajectory result. Another, more fundamental, limitation is that the length of the time window must be 
long enough to preserve an acceptable frequency resolution, even with the reassigned spectrogram. This limits the 
methods to weakly non-stationary signals. 

For signals with a rapid frequency modulation, the AML spectral estimation is well suited, as long as the noise is 
near iid. The size of the time window can be decreased because of the parametric nature of the method, since a priori 
knowledge has been taken into account. The performance is further increased by the use of a Kalman-like filter. The 
drawback is the necessity to find a good dynamical model for the evolution of the spectral components. We have 
chosen here the simplest model which works well for our experiments but the approach can be refined by increasing 
the number of parameters in order to consider more precisely the variation of the frequency. The AML algorithm also 
provides a quantitative estimation of the quality of the demodulation and the instantaneous power of the spectral 
component. Finally, the AML method has the advantage to provide directly frequency modulation as a function of 
time, in one stage. 
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AML ALGORITHM 

• First step : calculate Rj; using the following expressionil : 

= 20 E {x{^)X{^f+X{^)X{^f) , (32) 



with 



2Q 



X{i) = [x{i + K-l),x{i + K-2),..., x{i)]*'^ , (33) 



where * stands for complex conjugate. X is the complex conjugate of the time reversed version of X. 

• Second step : diagonalize R^^; one obtains the eigenvectors (Vi)i=i,,K and eigenvalues {Xi)i=i..K sorted in 
decreasing order. 

• Third step : Compute n„ and a^, using the set of equations 

K 

fln^ J2 ^^^^ ^ (34) 

i=M+l 

- ^T.(n„R.) ^ ^ E A.. (35) 

i=M+l 

• Forth step : choose F = F{t) as candidate value. Compute grad and H using 

grad= ^Re{Diag(S'+(F).n„(F).n„.S(F).P)} , (36) 
H = ^Rc {Diag ((S'+(F).n„(F).ri„.n„(F).S')) *P*} , (37) 
where the operator -k stand for the term to term matrix multiplication, and P* is the conjugate of P, and 

S' = [^,...,f^F. (38) 
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• Fifth step : using equations ( p7|) to (|30|), compute F(t + 1) and T(t + 1). 

The initiahzation of the algorithm is done by either (i) setting an initial value of F{1) or (ii) estimating this value 
using the maxima of the amplitude of the FFT of a small window of signal (of length 64 or 128 samples) and using 



the iterative algorithm described in section [II B 3 to converge towards F{1). 

For example, the extracted velocity of figure |^ is obtained by starting at the maximum of energy of the signal and 
applying the algorithm forward and backward in time. The algorithm is stopped as the mean of the inverse square 
root of the hessian over a window of size 400 samples exceeds 0.5 for more than 400 samples. 
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This is generally not the case, but this assumption remains valid as long as the loglikelihood is well approximated by its 
second order expansion around F. 



